Spreading of sexually transmitted diseases in heterosexual populations 
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The spread of sexually transmitted diseases (e.g. Chlamydia, Syphilis, Gonorrhea, HIV) across populations 
is a major concern for scientists and health agencies. In this context, both data collection on sexual contact net- 
works and the modeling of disease spreading, are intensively contributing to the search for effective immuniza- 
tion policies. Here, the spreading of sexually transmitted diseases on bipartite scale-free graphs, representing 
heterosexual contact networks, is considered. We analytically derive the expression for the epidemic threshold 
and its dependence with the system size in finite populations. We show that the epidemic outbreak in bipartite 
populations, with number of sexual partners distributed as in empirical observations from national sex surveys, 
takes place for larger spreading rates than for the case in which the bipartite nature of the network is not taken 
into account. Numerical simulations confirm the validity of the theoretical results. Our findings indicate that the 
restriction to crossed infections between the two classes of individuals (males and females) has to be taken into 
account in the design of efficient immunization strategies for sexually transmitted diseases. 



Disease spreading has been the subject of intense research 
since long time ago 1110, [H. On the one hand, epidemiolo- 
gists have developed mathematical models that can be used as 
a guide to understanding how an epidemic spreads and to de- 
sign immunization and vaccination policies IHI^llt]. On the 
other hand, data collections have provided information on the 
local patterns of relationships in a population. In particular, 
persons who may have come into contact with an infectious 
individual are identified and diagnosed, making it possible to 
contact-trace the way the epidemic spreads, and to validate 
the mathematical models. However, up to a few years ago, 
some of the assumptions at the basis of the theoretical mod- 
els were difficult to test. This is the case, for instance, of the 
complete network of contacts -the backbone through which 
the diseases are transmitted. With the advent of modern so- 
ciety, fast transportation systems have changed human habits, 
and some diseases that just a few years ago would have pro- 
duced local outbreaks, are nowadays a global threat for public 
health systems. A recent example is given by the severe acute 
respiratory syndrome (SARS), that spread very fast from Asia 
to North America a few years ago 0, HI 0]. Therefore, it is 
of utmost importance to carefully take into account as much 
details as possible of the structural properties of the network 
on which the infection dynamics occurs. 

Strikingly, a large number of statistical properties have been 
found to be common in the topology of real-world social, 
biological and technological networks 10, il- Of partic- 
ular relevance because of its ubiquity in nature, is the class 
of complex networks referred to as scale-free (SF) networks. 
In SF networks, the number of contacts or connections of a 
node with other nodes in the system, the degree (or connec- 
tivity) fc, follows a power law distribution, Pk ^ k L Recent 
studies have shown the importance of the SF topology on the 
dynamics and function of the system under study 13, 0]- 
For instance, SF networks are very robust to random failures, 
but at the same time extremely f rag ile to targeted attacks of 
the highly connected nodes jlOl [Till . In the context of dis- 
ease spreading, SF contact networks lead to a vanishing epi- 
demic threshold in the limit of infinite population when 7 < 3 



din [MB- 

This is because the exponent 7 is directly 
related to the first and second moment of the degree distri- 
bution, (fc) and (fc^), and the ratio (fc)/(fc^) determines the 
epidemic threshold above which the outbreak occurs. When 
2 < 7 < 3, (fc) is finite while (fc^) goes to infinity, that 
is, the transmission probability required for the infection to 
spread goes to zero. Conversely, when 7 > 3, there is a finite 
threshold and the epidemic survives only when the spreading 
rate is above a certain critical value. The concept of a criti- 
cal epidemic threshold is central in epidemiology. Its absence 
in SF networks with 2 < 7 < 3 has a number of impor- 
tant implications in terms of prevention policies: if diseases 
can spread and persist even in the case of vanishingly small 
transmission probabilities, then prevention campaigns where 
individuals are randomly chosen for vaccination are not much 
effective |[Il[Il[ll[ll. 

Our knowledge of the mechanisms involved in disease 
spreading as well as on the relation between the network struc- 
ture and the dynamical patterns of t he sp readin g proc ess has 
improved in the last several years III6I [TtI M8L M9D. Cur- 



rent approaches are either individual-based simulations 111811 
or metapopulation models where network simulations are car- 
ried out through a detailed stratification of the population 
and infection dynamics ll20ll . In the particular case of sexu- 
ally transmitted diseases (STDs), infections occur within the 
unique context of sexual encounters, and the network of con- 
tacts [[illllllllllllHlllllilisa critical ingredient of any 
theoretical framework. Unfortunately, ascertaining complete 
sexual contact networks in significatively large populations is 
extremely difficult. However, here we show that it is indeed 
possible to make use of known global statistical features to 
generate more accurate predictions of the critical epidemic 
threshold for STDs. 



I. NETWORKS OF SEXUAL CONTACTS 



Data from national sex surveys 112 IL l22l I23L l24l 12511 provide 
quantitative information on the number of sexual partners, the 
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degree fc, of an individual. Usually, surveys involve a random 
sample of the population stratified by age, economical and 
cultural level, occupation, marital status, etc. The respondents 
are asked to provide information on sexual attitudes such as 
the number of sex partners they have had in the last 12 months 
or in their entire life. Although in most cases the response rate 
is relatively small, the information gathered is statistically sig- 
nificant and global features of sexual contact patterns can be 
extracted. In particular, it turns out that the number of hetero- 
sexual partners reported from different populations is well de- 
scribed by power-law SF distributions. TableUsummarizes the 
main results of surveys conducted in Sweden, United King- 
dom, Zimbabwe and Burkina Faso ll2lll22i l23ll24ll . 

The first thing to notice is the gend er-sp ecific difference in 
the number of sexual acquaintances l2ll I22I l23l l24ll . This 
is manifested by the existence of two different exponents in 
the SF degree distributions, one for males (7^/) and one for 
females (7_f). Interestingly enough, the predominant case in 
Table m no matter whether data refers to time frames of 12 
months or to entire life, consists of one exponent being smaller 
and the other larger than 3. This is certainly a borderline case 
that requires further investigation on the value of the epidemic 
threshold. 

The differences found in the two exponents and 7^/ 
have a further implication for real data and mathematical mod- 
eling. In an exhaustive survey, able to reproduce the whole 
network of sexual contacts, the total number of female part- 
ners reported by men should equal the total number of male 
sexual partners reported by women. Mathematically, this 
means that the number of links ending at population M (of 
size Nm) equals the number of links ending at population F 
(of size Np), which translates into the following closure rela- 
tion; 

NF{k)F = NM{k)M ■ (1) 

Assuming that the degree distributions for the two sets are 
truly scale-free, then = • k^'^'^, with the symbol 

G standing for the gender (G = F, M), and fco being the 
minimum degree. Moreover, if Nq >> 1 and 70 > 2 for any 
G, Eq. ([T]i gives the relation between the two population sizes 
as 

= Nf§^ ^ nJ2^) i2) 

which implies that the less heterogeneous (in degree) popula- 
tion must be larger than the other one. 

In conclusion, the empirical observation of two different 
exponents demands for a more accurate description of the 
network of heterosexual contacts as bipartite SF graphs, i.e. 
graphs with two set of nodes, and links connecting nodes from 
different sets only. In the following we will consider a graph 
with Nm nodes, representing males, characterized by the ex- 
ponent 7j\f , and Nm nodes, representing females, character- 
ized by 7i?. Concerning the choice of the couple of exponents 
from those reported in Table|I] one must be careful that differ- 
ent STDs have different associated (recovery) time scales, and 
that the spreading is based on the assumption that the links are 




FIG. 1 : Bipartite and unipartite networks.(a) Representation of a 
bipartite network accounting for heterosexual contact networks. In 
such a network we have Nm (Nf) nodes representing males (fe- 
males), and only male-female links are allowed. Figure (b) repre- 
sents a rewired version of network (a) where the bipartite nature is 
lost, while the degree of the nodes is preserved. The two graphs have 
the same couple of degree distributions (one for males and one for 
females), although only (a) reflects the bipartite character of hetero- 
sexual contact networks. 



concurrent on the time scale of the disease. In this sense, the 
exponents extracted from one-year data seem better suited to 
most of the STDs, being HIV an important exception. How- 
ever, during the lifetime of sexually active individuals, sex- 
ual behavior is likely to change due to changes in residence, 
marital status, age-linked sexual attitudes, etc We thus 
prefer to use life-cycle data collections, that integrate all these 
patterns and can consequently be regarded as better statistical 
indicators. After all, the values reported in Table U indicate 
that both one-year and cumulative data produce exponents in 
the same range. 



II. THEORETICAL MODELING 

The problem of how a disease spreads in a population con- 
sisting of two classes of individuals can be tackled by invoking 
the so-called criss-cross epidemiological model fl]. As illus- 
trated in the bipartite network of Fig.[T](a), in the criss-cross 
model, the two populations of individuals (N^^ males and N^ 
females) interact so that the infection can only pass from one 
population to the other by crossed encounters between the in- 
dividuals of the two populations, incorporating in this way 
one of the basic elements of the heterosexual spreading of 
STDs ll28ll . In particular, we consider a susceptible-infected- 
susceptible (SIS) dynamics, in which individuals can be in 
one of two different states, namely, susceptible (S) and infec- 
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TABLE I: Statistical properties of sexual contact networks from national sex surveys conducted in four different countries: Sweden, United 
Kingdom, Zimbabwe and Burkina Faso. The exponents 7f and 7^/ are referred to the distribution of number of sexual partners cumulated in 
12 months and in the respondent's lifetime. The number of respondents is also reported. 



survey ref. 


7f (12 months) 7A/ (12 months) 


7F (life) 7A/ (life) 


Resp. Resp. (F) Resp. (M) 


Sweden [21] 


3.54 ±0.20 3.31 ±0.20 


3.1 ±0.30 2.6 ±0.30 


2810 


U.K. 122,23} 


3.10 ±0.08 2.48 ±0.05 


3.09 ±0.20 2.46 ±0.10 


11161 6399 4762 


Zimbabwe [231 


2.51 ±0.40 3.07 ±0.20 


2.48 ±0.15 2.67 ±0.18 


9843 5424 4419 


Burkina Faso [241 


3.9 ±0.2 2.9 ±0.1 




466 179 287 



tious (X). If S'^' and l'^^ (S^, I^) stand for a male (female) 
respectively in the susceptible and in the infectious state, the 
epidemic in the SIS criss-cross model propagates by the fol- 
lowing mechanisms: 



5^- 



jF '^M^ jM 

^ S^, 



j-M MA/ 



s 



M 



being vm, vp, and /^ip the infection and recovery prob- 
abilities for males and females. In the case of heterogeneous 
contact networks, there is a further compartmentalization of 
the population into classes of individuals with the same de- 
gree k, i.e. the same number of sexual partners. Denoting the 
fraction of males (females) with degree k in the susceptible 



or infectious state by sjr and ij^ (si and i^j, respectively, 
and adopting a mean field approach |l3l[i2|,[i3lj the differen- 
tial equations describing the time evolution of the densities of 
susceptible and infected individuals in each population are 



1 di^it) 



dt 



-if(t)±A^^fc[l-zf(t)]ef(t), (3) 



^J.M 



dt 



M 



k [1 



(4) 



where, \g ~ vqI [iq (G = M) are the effective transmis- 
sion probabilities. The quantities Q^^t) , Q^{t) stand for the 
probabilities that a susceptible node of degree k of one pop- 
ulation encounters an infectious individual of the other set. 
Equations (|3]l and (HJlhave the same functional form of the 
equation derived in ||12I1 for unipartite networks. Neglecting 
degree-degree correlations, the critical condition for the oc- 
currence of an endemic state reduces to: 



V A_fAa/ > a* — 



{k)F{k}M 



(5) 



fact, working with a unipartite representation of a sexual net- 
work, as that shown in Fig.[T](b), with Nm + Np nodes and a 
degree distribution = {NmP^ + NpP[)/{Nm + Np), 
one can express the epidemic threshold as a function of the 
first and second moments of the male and female degree dis- 
tributions as 



(fc2) {k^)Al{k)F + {k^)p{k)M 



2{k)M{k)l 



(6) 



Equations ^ and ^ are clearly different. For real SF net- 
works of sexual contacts, the two thresholds are finite (in the 
infinity size limit) only when the two exponents 7^/ and 'yp 
are both larger than 3, e.g. for the one-year number of partners 
in Sweden (see Table IJi. In such a case, the two expressions 
read: 



K 

Xr 



l_ / (7F-3)(7A/-3) 
fcoV (7F-2)(7M-2) 



(7) 



2(7F-3)(7M-3) 

ko [(7F - 2)(7M - 3) ± (7A/ - 2)(7F - 3)| 



These two thresholds are only equal in the case 7m = 7f 
studied in ifTsll . More importantly, when jp =^ 7m, we have 
X'c > Ac, that is the epidemic state in bipartite networks oc- 
curs for larger transmission probabilities than in unipartite 
networks. This result is good news and highlights the im- 
portance of incorporating the crossed infections scheme in the 
propagation of STDs. However, as shown in Table U most of 
the real networks have at least one exponent 7^ (G = F, M) 
smaller than 3. This means that, in most of the practical cases, 
the two epidemic thresholds vanish as the system size goes to 
infinity, no matter the formulation used to model the disease 
propagation. On the other hand, real populations are finite and 
thus the degree distributions have a finite variance regardless 
of the exponents. Consequently, an epidemic threshold does 
always exist and, in order to compare unipartite with bipartite 
networks, one must then pay attention to the scaling of the 
threshold with the size of the population. 



yielding that a necessary condition for the absence of the epi- 
demic threshold is the divergence of at least one of the second 
moments of the degree distributions, (fc^)M and {k"^) p. Equa- 
tion ^ can be compared with the condition obtained without 
taking into account that, in heterosexual networks, the infec- 
tion can occur only between male-female couples lfT2i [l3ll . In 



A. Finite Populations 

We now analyze in more details the differences between Ac 
and Ac when one of the two exponents (let us say 7^/ without 
loss of generality) is in the range 2 < 7m < 3, while ^p > 2. 
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FIG. 2: Epidemic tliresliolds as a function of tlie total population size. The thresholds A* (Eq. ^) and Ac (Eq. {SJ) obtained for bipartite 
and unipartite networks, are plotted as a function of the population size A'' for two networks with degree distributions as those found for the 
sexual networks of (a) Sweden (2lll and (b) United Kingdom |(23l. The two insets show that the ratio A*/Ac grows as increases, so that for 
a typical population size of A*" = lO'' the thresholds for heterosexual networks are respectively 53% (Sweden) and 130% (U.K.) larger than 
the values expected for the unipartite networks. 



First, we derive the size scaling of the critical threshold in 
unipartite graphs, Ac. Equation ^ yields 

Np{k)F + NM{k)M 



Np(k^)p + NF{k^)M ' 

Manipulating this expression by considering again the limit of 
large (but finite) population sizes, Nq >> 1 {G = F, A/), we 
obtain 

, ^, 2(3 - 7m)/(7m - 2) 



fco 



3-7M -IF -2 
3 — 1F 1M—2 



The final expression for Ac can now be obtained by using the 
closure relation of Eq. (|2]i for the M and F population sizes, 
yielding: 

X ^, 2(3 - 7m)/(7m - 2) 



In this formula, only one population size Nm appears. Finally, 
if e.g. > lAi, the above equation reduces to 



2(3 "7m) ^^ 



fco(7M - 2) 



(9) 



that contains simultaneously the cases when 2 < 7^ < 3 and 

7f > 3. 

Now we calculate the scaling of the epidemic threshold in 
bipartite (heterosexual) networks. A*. Manipulating Eq. (|5]), 
A* can be as well expressed as a function of the two exponents 



TABLE II: Scaling exponents of the epidemic thresholds.Scaling 

exponents, a and a*, of the epidemic thresholds, Ac ~ A'^^ ^rid 
Ac ~ N^j , obtained for the SIS model on unipartite networks and 
when a bipartite network is considered, respectively. The two sit- 
uations considered (2 < 7f < 3 and 7f > 3) conespond to 
2 < 7M < 3. 



Network 


a* 


a 


2 < 7F < 3 




f 3 IF 1 3-7Af ^ 




3-7i\J 
7J\/-1 


7F > 3 


1 


( S-IM ] 




3 — IM 
7J\/-1 



7A/ and 7^-, and one population size: 



B 



1 



7F-2 7M-1 ^ 
7j,/~2 jF-l " 



with B = , '^'"'^ \ ■ The above expression, when eval- 

'=o(2-7m)(2-7f) ^ 

uated for 2 < 7^ < 3 (G = i^, Af ) and, e.g., •jf > Jm yields 



B 



1/2 /tf^Tm^A^^ Kif^H-^^ 



7a/ - 2 7f - 1 



M 



(10) 

On the other hand, when e.g. jf > 3, the expression reduces 
to 



' (3-7m)(7f^3) 
(2-7A/)(2-7F)fco' 



(11) 
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FIG. 3: Monte Carlo simulations for the UK and Sweden lifetime networks.The SIS phase diagram, (/(t)) vs. A/A*, is reported for 
synthetic networks of different sizes. Note that for each curve the x-axis has been rescaled by the theoretical value (Eq. [S} for the critical point 
in bipartite networks. A*, of the corresponding network size and exponents 7f and jm- These exponents used in the network construction (see 
Methods) are those extracted from the lifetime degree distributions (see Table|l]l of Sweden I21I1 . panels (a) and (c), and U.K. (23.l23ll . panels 
(b) and (d). The results plotted in (a) and (b) correspond to the SIS model on a bipartite network, whereas those shown in (c) and (d) correspond 
to a unipartite substrate. The numerical results clearly indicate (see insets) the validity of the analytical predictions for the epidemic threshold 
in heterosexual (bipartite) contact networks and the underestimation of the epidemic threshold when a unipartite substrate is employed. 



B. Comparing the scalings 



III. NUMERICAL SIMULATIONS 



Although both epidemic thresholds, Ac and A*, tend to 
zero as the population goes to infinity, the scaling relations, 

XciNM,jM,jF) Nj^r and A*(^m, 7a/, 7i^) - A^m"'- 
are characterized by two different exponents, a and a*. Table 
|II]reports the expression of these two exponents as a function 
of and "fM, showing that a* is always smaller than a. In 
particular, for the most common case (see Table i.e. when 
one degree distribution exponent is in the range ]2, 3], and the 
other one is larger than 3, the value of a* found for bipartite 
networks is two times smaller than a. As a consequence, the 
results show that in finite bipartite populations the onset of the 
epidemic takes place at larger values of the spreading rate. In 
other words, it could be the case that for a given transmission 
probability, in the unipartite representation shown in Fig.[T](b) 
the epidemic would have survived infecting a fraction of the 
population, while when only crossed infections are allowed, 
as in Fig.[T](a), the same disease would not have produced an 
endemic state. Moreover, the difference between the epidemic 
thresholds predicted by the two approaches increases with the 
system size. This dependency is shown in Fig. |2] where we 
have reported, as a function of the system size, the critical 
thresholds obtained by numerically solving Eqs. ^ and (|6]l 
with the values of and -fp found for the lifetime distribu- 
tion of sexual partners in Sweden ll2lll and U.K ll22[|23ll . 



To check the validity of the analytical arguments and also 
to explore the dynamics of the disease above the epidemic 
threshold, we have conducted extensive numerical simula- 
tions of the SIS model in bipartite and unipartite computer- 
generated networks. Bipartite and unipartite graphs of a given 
size are built up (see Methods section) having the same de- 
gree distributions, and P^, and thus they only differ in 
the way the nodes are linked. A fraction of infected individ- 
uals is initially randomly placed on the network and the SIS 
dynamics is evolved: at each time step susceptible individuals 
get infected with probability u if they are connected to an in- 
fectious one, and get recovered with probability /x = 1 (hence, 
the effective transmission probability is A = v). After a tran- 
sient time, the system reaches a stationary state where the total 
prevalence of the disease, {I{t)), is measured (see Methods). 
The results are finally averaged over different initial condi- 
tions and network realizations. Fig. |3] shows the fraction of 
infected individuals as a function of A/ A* for several system 
sizes and for the bipartite ((a) and (b)) and unipartite ((c) and 
(d)) graphs. In this figure, the infection probability A has been 
rescaled by the theoretical value A* given by Eq. The pur- 
pose of the rescaling is twofold. First it allows to check the 
validity of the theoretical predictions and, at the same time, 
it provides a clear comparison of the results obtained for bi- 
partite networks with those obtained for the unipartite case. 
Again we have used the values of jm and jp extracted from 
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the lifetime number of sexual partners reported for Sweden 
and U.K. iUllllll. Fig.[3]indicates that the analytical solu- 
tion, Eq. dD, is in good agreement with the simulation results 
for the two-gender model formulation. Conversely when the 
bipartite nature of the underlying graph is not taken into ac- 
count, the epidemic threshold is underestimated, being Ac/A* 
smaller than 1. In addition to this, the error in the estimation 
grows as the population size increases, in agreement with our 
theoretical predictions. 

IV. CONCLUSIONS 

The inclusion of the bipartite nature of contact networks to 
describe crossed infections in the spread of STDs in hetero- 
sexual populations is seen to affect strongly the epidemic out- 
break and leads to an increase of the epidemic threshold. Our 
results show that, even in the cases when the epidemic thresh- 
old vanishes in the infinite network size limit, the epidemic 
incidence in finite populations is less dramatic than actually 
expected for unipartite scale-free networks. The results also 
point out that the larger the population, the greater the gap be- 
tween the epidemic thresholds predicted by the two models, 
therefore highlighting the need to accurately take into account 
all the available information on how heterosexual contact net- 
works look like. Our results also have important consequences 
for the design and refinement of efficient degree-based im- 
munization strategies aimed at reducing the spread of STDs. 
In particular, they pose new questions on how such strategies 
have to be modified when the interactions are further compart- 
mentalized by gender and only crossed infections are allowed. 
We finally stress that the present approach is generalizable 
to other models for disease spreading {e.g. the ''susceptible- 
infected-removed" model) and other processes where crossed 
infection in bipartite networks is the mechanism at work. 

V. METHODS 
A. Bipartite Network construction 

Synthetic bipartite networks construction starts by fixing 
the number of males, Nm and the two exponents 7^/ and 
of the power-law degree distributions corresponding to males 
and females respectively. The first stage consists of assign- 
ing the connectivity kf^ (i = 1,...,Nm) to each member of 
the male population by generating random numbers with 
probability distribution P^^ = ^7,//?"^^^ (J2To ^A/fc"'^" = 
1, with fco = 3). The sum of these N„i random numbers 
fixes the number of links Ni of the network. The next step 
is to construct the female population by means of an iterative 
process. For this purpose we progressively add female indi- 



viduals with a randomly assigned degree following the distri- 
bution Pf = Apk^'^'' Apk-'^^ = 1, with fco = 3). 
Female nodes are incorporated until the total female connec- 
tivity reaches the number of male edges, kf < iV;. In this 
way one sets the total number of females Np- Once the two 
sets of Nm males and Np females with their corresponding 
connectivities are constructed each one of the Ni male edges 
is randomly linked to one of the available female edges avoid- 
ing multiple connections. Finally those few female edges that 
did not receive a male link in the last stage are removed and 
the connectedness of the resulting network is checked. 

B. Unipartite Network construction 

Synthetic unipartite networks has been constructed in two 
ways. The simplest one consists of taking the two sets of Nm 
males and Np females constructed for the bipartite network 
and apply a rewiring process to the entire population, i.e. al- 
lowing links between individuals of the same sex. In the sec- 
ond method, a set of iV = Nm + Np individuals whose con- 
nectivities are randomly assigned following the degree dis- 
tribution P{k) = {Nm/N)P^^ + {Np/N)P[ is generated 
before applying a wiring process between all pairs of edges. 
In both methods the wiring process avoids multiple and self 
connections and those isolated edges that remain at the end of 
the network construction are removed. The connectedness of 
the networks is also checked. 



C. Numerical Simulations of SIS dynamics 

Montecarlo simulations of SIS dynamics are performed us- 
ing networks of sizes ranging from = 2 • 10"* to iV = 8-10*. 
The initial fraction of infected nodes is set to 1% of the net- 
work size. The SIS dynamics is initially evolved for a time 
typically of 10''^ time-steps and after this transient the system 
is further evolved over consecutive time windows of 2 • 10'^ 
steps. In these time windows we monitor the mean value of 
the number of infected individuals, (/(t)). The steady state is 
reached if the absolute difference between the average num- 
ber of infected individuals of two consecutive time windows 
is less than 1/^/N. 
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